1.1 Packages

Load required packages into R session.

library(devtools)
library(tidyverse)
library(Biostrings)

1.2 Functions

url <- "https://raw.githubusercontent.com/hickslab/ProgenesisLFQ/master/"

source_url(paste0(url, "R/ProgenesisLFQ_Global.R"))

## SHA-1 hash of file is 46d72ef9fbf7843a06a8c125d64eaf3c02bfa2e5

1.3 Workflow

Load Protein Measurements spreadsheet exported from Progenesis.

protm <- "data/20180502_WOS52_Cr_UPS_protm.csv" %>%
  paste0(url, .) %>%
  read_csv(., skip = 2, col_types = cols())

Separte leading protein accession from group members.

• Some protein accessions were grouped together by semicolons (;) in the data. These groups represent homologous proteins with shared peptide identifications. The first accession in each group is generally the most confident inference and is used for analysis, other accessions are considered group members.

protm <- protm %>%
  separate(Accession,
           into = c("Accession", "Group members"),
           sep = ";",
           extra = "merge",
           fill = "right")

protm

Define column indeces for the normalized abundance values.

• The Protein Measurements export from Progenesis contains the normalized and raw abundances for each raw file in the experiment.

protm %>% names()

##  [1] "Accession"              "Group members"         
##  [3] "Peptide count"          "Unique peptides"       
##  [5] "Confidence score"       "Anova (p)"             
##  [7] "Max fold change"        "Highest mean condition"
##  [9] "Lowest mean condition"  "Description"           
## [11] "20141222_WOS521"        "20141222_WOS526"       
## [13] "20141222_WOS5211"       "20141222_WOS5216"      
## [15] "20141222_WOS522"        "20141222_WOS527"       
## [17] "20141222_WOS5212"       "20141222_WOS5217"      
## [19] "20141222_WOS523"        "20141222_WOS528"       
## [21] "20141222_WOS5213"       "20141222_WOS5218"      
## [23] "20141222_WOS521_1"      "20141222_WOS526_1"     
## [25] "20141222_WOS5211_1"     "20141222_WOS5216_1"    
## [27] "20141222_WOS522_1"      "20141222_WOS527_1"     
## [29] "20141222_WOS5212_1"     "20141222_WOS5217_1"    
## [31] "20141222_WOS523_1"      "20141222_WOS528_1"     
## [33] "20141222_WOS5213_1"     "20141222_WOS5218_1"

samples <- 11:22

samples

##  [1] 11 12 13 14 15 16 17 18 19 20 21 22

Filter to remove proteins from the contaminant database.

Filter to remove proteins if not enough peptide evidence.

• Proteins are identified by unique and shared PSMs, so usually the confidence in our protein identification is greatest with increasing unique peptides. It is common in proteomics to remove “one-hit-wonders” by requiring > 1 unique peptide per protein.

Select the identifier and abundance columns.

data <- protm %>%
  filter(Description != "cRAP") %>%
  filter(`Peptide count` >= 2 & `Unique peptides` >= 1) %>%
  select(Accession, samples) %>%
  data.frame()

data

Example plot construction.

plot_jitter <- function(df){
  temp.df <- df %>%
    gather(sample, abundance, -1)
  
  temp.df %>%  
    ggplot(., aes(x = sample, y = abundance, color = sample)) +
    geom_jitter(alpha = 0.5) +
    geom_boxplot(color = "black",
                 fill = NA,
                 outlier.shape = NA,
                 size = 1.5) +
    guides(color = FALSE, fill = FALSE) +
    theme_classic(base_size = 32) +
    coord_flip()

}

data %>% plot_jitter()

2.1 Packages

Load required packages into R session.

library(devtools)
library(tidyverse)
library(Biostrings)

2.2 Functions

url <- "https://raw.githubusercontent.com/hickslab/ProgenesisLFQ/master/"

source_url(paste0(url, "R/ProgenesisLFQ_Peptide.R"))

## SHA-1 hash of file is cc86fbeb8fda5d978f0e08b1448e42efc52d730d

2.3 Workflow

pepm <- "data/20180502_WOS52_Cr_UPS_pepm.csv" %>%
  paste0(url, .) %>%
  read_csv(., skip = 2, col_types = cols())

protm <- "data/20180502_WOS52_Cr_UPS_protm.csv" %>%
  paste0(url, .) %>%
  read_csv(., skip = 2, col_types = cols())

pepm <- "data/20190123_EWM_AZD1_R_rank-lessthan11-include_uniprot_pepm.csv" %>%
  paste0(url, .) %>%
  read_csv(., skip = 2, col_types = cols())

pepm

protm <- "data/20190123_EWM_AZD1_R_rank-lessthan11-include_uniprot_protm.csv" %>%
  paste0(url, .) %>%
  read_csv(., skip = 2, col_types = cols()) %>%
  select(1:4) %>%
  separate_rows(., Accession, sep = ";")

protm

database <- "data/Cr_uniprot_crap_20190130.fasta" %>%
  paste0(url, .) %>%
  readAAStringSet(.)

names(database) <- str_split(names(database), " ", simplify = TRUE)[, 1]

database

##   A AAStringSet instance of length 18944
##         width seq                                      names               
##     [1]   751 MTISTPEREAKKVKIAVDR...GGIATTWSFFLARIISVG sp|P12154|PSAA_CHLRE
##     [2]   735 MATKLFPKFSQGLAQDPTT...YIFTYAAFLIASTSGRFG sp|P09144|PSAB_CHLRE
##     [3]    81 MAHIVKIYDTCIGCTQCVR...SVRVYLGSESTRSMGLSY sp|Q00914|PSAC_CHLRE
##     [4]    43 MIFDFNYIHIFMLTITSYV...LVFTLGIYLGLLKVVKLI sp|P50369|PETL_CHLRE
##     [5]   160 MSVTKKPDLSDPVLKAKLA...LGIGSTFPIDISLTLGLF sp|P23230|PETD_CHLRE
##     ...   ... ...
## [18940]   238 MSKGEELFTGVVPILVELD...LEFVTAAGITHGMDELYK sp|GFP_AEQVI
## [18941]   204 MAEEVEEERLKYLDFVRAA...LPLLPTEKITKVFGDEAS sp|SRPP_HEVBR
## [18942]   138 MAEDEDNQQGQGEGLKYLG...SSLPGQTKILAKVFYGEN sp|REF_HEVBR
## [18943]   348 MFSSVMVALVSLAVAVSAN...VMNADNHEYFSENNPAQS sp|PLMP_GRIFR
## [18944]   271 MSHIQRETSCSRPRLNSNL...DNPDMNKLQFHLMLDEFF sp|KKA1_ECOLX

#pepm %>% names()

samples <- 19:34

samples

##  [1] 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

data <- pepm %>%
  filter(Score > 13) %>%
  filter(Description != "cRAP") %>%
  left_join(., protm, by = "Accession") %>%
  #reduce_features() %>%
  #filter_redox(., reduced = "NEM") %>%
  reduce_features_tidy() %>%
  filter_redox_tidy(reduced = "Nethylmaleimide") %>%
  get_identifier_redox(database = database, reduced = "NEM") %>%
  reduce_identifiers(group = samples) %>%
  select(Identifier, samples) %>%
  data.frame()

## Warning: slice_() is deprecated. 
## Please use slice() instead
## 
## The 'programming' vignette or the tidyeval book can help you
## to program with slice() : https://tidyeval.tidyverse.org
## This warning is displayed once per session.

## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.

data

3.1 Packages

Load required packages into R session.

library(imp4p)
library(broom)

3.2 Functions

url <- "https://raw.githubusercontent.com/hickslab/ProgenesisLFQ/master/"

source_url(paste0(url, "R/StatLFQ.R"))

## SHA-1 hash of file is 21bea41dbe4f674300b2d84ef64282ad79962cf0

3.3 Workflow

Define input data.

data <- "data/20180502_WOS52_Cr_UPS_protm_Process.csv" %>%
  paste0(url, .) %>%
  read_csv(., col_types = cols()) %>%
  data.frame()

Define the column indeces for replicates in each condition.

• Assuming input data has simplified format of an identifier column followed by abundance columns for each raw file in the experiment.

a <- 2:5; b <- 6:9; c <- 10:13 # ???

group <- list("25" = a, "50" = b, "100" = c) # ???

group.compare <- list("25-50" = list(a, b),
                      "25-100" = list(a, c),
                      "50-100" = list(b, c))

a

## [1] 2 3 4 5

b

## [1] 6 7 8 9

c

## [1] 10 11 12 13

group

## $`25`
## [1] 2 3 4 5
## 
## $`50`
## [1] 6 7 8 9
## 
## $`100`
## [1] 10 11 12 13

group.compare

## $`25-50`
## $`25-50`[[1]]
## [1] 2 3 4 5
## 
## $`25-50`[[2]]
## [1] 6 7 8 9
## 
## 
## $`25-100`
## $`25-100`[[1]]
## [1] 2 3 4 5
## 
## $`25-100`[[2]]
## [1] 10 11 12 13
## 
## 
## $`50-100`
## $`50-100`[[1]]
## [1] 6 7 8 9
## 
## $`50-100`[[2]]
## [1] 10 11 12 13

Rename the abundance columns in a simplified “Condition-Replicate” format.

data <- data %>%
  rename_columns(., group)

data %>% names()

##  [1] "Accession" "25-1"      "25-2"      "25-3"      "25-4"     
##  [6] "50-1"      "50-2"      "50-3"      "50-4"      "100-1"    
## [11] "100-2"     "100-3"     "100-4"

Filter to keep identifiers with >50% of replicates having nonzero abundances in any condition.

Perform a log₂-transformation of abundances as a variance-stabilization. The base R function log2() returns “-Inf” for zero values. These instances are changed to “NA” following the transformation.

Impute missing values (“NA”) using a conditional strategy with the imp4p package. Iterate by condition and check if each idenfifier has reliable quantitation.

• If at least one replicate has a nonzero abundance, impute other replicates with values drawn from a normal distribution centered on the mean of nonzero replicates.

• If all replicates have nonzero abundance, impute with small values drawn from a normal distribution centered on the lower 25^th percentile of abundances.

data2 <- data %>%
  clean_min(., group, nonzero = 3) %>%
  transform_data(., group, method = "log2") %>%
  impute_imp4p(., group)

data2 %>% plot_jitter()

data3 <- data2 %>%
  calculate_ttest(., group.compare, fdr = TRUE)

data4 <- data2 %>%
  calculate_1anova()

data3 <- data3 %>%
  calculate_fc(., group.compare, difference = TRUE)


data4 <- data4 %>%
  calculate_fc(., group.compare, difference = TRUE)

#data4 %>% calculate_hclust(., group, k = 3)

4.1 Packages

Load required packages into R session.

4.2 Functions

url <- "https://raw.githubusercontent.com/hickslab/ProgenesisLFQ/master/"

source_url(paste0(url, "R/PlotLFQ.R"))

## SHA-1 hash of file is 3210231e1ee84ce77a0cbb73f2ed57b20bfa6d93

4.3 Theme

Define custom plotting elements to create a theme.

theme_custom <- function(base_size = 32){
  theme_bw(base_size = base_size) %+replace%
    theme(
      strip.background = element_blank(),
      #strip.text.x = element_text(size = 18),
      #strip.text.y = element_text(size = 18),
      #axis.text.x = element_text(size=14),
      #axis.text.y = element_text(size=14,hjust=1),
      axis.ticks =  element_line(colour = "black"), 
      #axis.title.x= element_text(size=16),
      #axis.title.y= element_text(size=16,angle=90),
      panel.background = element_blank(), 
      panel.border = element_blank(), 
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(), 
      #panel.margin = unit(1.0, "lines"), 
      plot.background = element_blank(), 
      plot.margin = unit(c(0.5,  0.5, 0.5, 0.5), "lines"),
      axis.line.x = element_line(color="black", size = 1),
      axis.line.y = element_line(color="black", size = 1)
    )
}

theme_custom2 <- function(){
  theme_classic(base_size = 32) %+replace%
    theme(
      strip.background = element_blank()
      
    )
}

4.4 PCA

Principal component analysis (PCA)

data2 %>% plot_pca(., group)# + theme_custom()

4.5 Volcano Plot

Biological versus statistical significance

data3 %>%
    plot_volcano(.,
                 group,
                 group.compare,
                 fdr = TRUE,
                 threshold = 2,
                 xlimit = 8,
                 ylimit = 8) +
  theme_custom()

4.6 Trend Profiles

Visualization of results from heirarchical clustering.

data4 %>%
  filter(fdr < 0.05) %>%
  #filter(abs(`0-60_FC`) >= 1) %>%
  plot_hclust(., group, k = 2) +
  theme_custom()

5.1 Packages

Load required packages into R session.

5.2 Functions

url <- "https://raw.githubusercontent.com/hickslab/ProgenesisLFQ/master/"

source_url(paste0(url, "R/AnnotateLFQ.R"))

## SHA-1 hash of file is 34d87fddd0d3e4a15a837ac42bbc820b3df5a648

5.3 UniProtKB API

Given a list of protein accessions, access UniProtKB and pull known information.

6.1 Dataframes

The write_csv() function from the readr package in tidyverse is one option to save a dataframe copy in the working directory.

#data5 %>% write_csv(., "data5.csv")

6.2 Plots

Consistent plot size/resolution.

#pdf("figure.pdf", width = 10, height = 10)
#png("figure.png", width = 12, height = 9, units = "in", res = 300)
# RUN PLOT
#dev.off()